＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
■■■■
■ ■ ■ ■
■ ■■■ ■■■
■ ■ ■ ■
■■■■
〜基礎から ★ Ｃ＋＋Ｐｒｏｇｒａｍｉｎｇ〜
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
【注意】 このマガジンは、最大化してお読みください。
また、等角フォントでお読みください。
（ＭＳ ゴシックなど）
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
発行者 むーくん
マガジンＮＯ． アルゴリズム編No.54 ヤコビ法
発行日 02/10/04
講読人数 ３８００名ぐらい
マガジンＩＤ 0000050494
このマガジンは、まぐまぐから配信されています。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★あいさつ★
お久しぶりです ^^;
なんだか、No.53が発送されてなかったみたいで、
随分間が開いてしまいました・・・
ああ、不手際でした・・・
今月も頑張りますのでこれからもよろしくお願いいたします。m(--)m
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
■退職金くじ１０万円【この場】で当る！ http://www.jobfront.to/a.html
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★目次★
・まえがき
・ヤコビ(Jacobi)法
・サンプルプログラム
・今日のポイント
・予告
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★まえがき★
固有値・固有ベクトルを求めるためには、固有方程式という
代数方程式をとき、その後掃き出し法を適用するという面倒な物でした。
ところが、対称となる行列が「対称行列」の場合、
もっと簡単な方法が存在します。
今回はその中の最も有名で簡単な方法、「ヤコビ法」を学習します。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★対称行列と固有行列★
対称行列とは、
| a b c |
| b e f | の様に、対角成分（aeiのライン）を境として対称になっている
| c f i | のことを指します。
対称行列の固有値、λ1,λ2,λ3・・・を対角成分に並べた
行列は常に成分が実数値をとり、
固有ベクトルx1,x2,x3・・・を縦に並べた行列は、
直交行列になるという性質があります。
直交行列の性質は以前学習したように、
T[P] = P^-1
と、転置行列と逆行列が等しくなります。
＃Pは、固有ベクトルを縦に並べたものです。
固有値の定義、AP = Pλの両辺に左からP^-1を掛けることによって
P^-1・A・P = λ が得られます。
λは固有値を対角成分に並べた対角行列です。
| λ1 0 0 |
＃つまり、 λ = | 0 λ2 0 | ということです。
| 0 0 λ3 |
普通は大文字で「Λ」と書くことが多いですが
お気になさらずに。
これを「対角化」というのでした。
ここで、対称行列の場合、固有ベクトルを縦に並べた行列が直交行列なので
T[P] = P^-1 が利用でき、
T[P]AP = λ と簡単化することができます。
ここまでは前回までの復習です。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★ヤコビ(Jacobi)法★
さて、上のように、「ピッタリ正しいPを使えば」、
λが求まり、対角化が可能な訳ですが、ここで逆を考えてみます。
適当な直交行列 Uを当てはめて、T[U]AUが対角行列なら、
Tは固有ベクトル、λは固有値と考えて差し支えないわけです。
直交行列は適当にとるので、一発で一致することは少ないですが、
何回も繰り返すことで真値に漸近させることができるというわけです。
最初に与えられた行列AをA(0)、m回直交変換した行列をA(m)と
書くことにします。すると、
t[Un]t[Un-1]・・・t[U2]t[U1]A[U1][U2]・・・[Un-1][Un] = A(n)
簡単に書けば
t[U](n+1)A(n)[U](n+1) = A(n+1)
と書けますね。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★直交行列の選び方★
実際に問題となってくるのは直交行列の選び方です。
以前学習したとおり、
| cosθ sinθ |
| -sinθ cosθ | は直交行列であると学習しました。
ヤコビ法では、以下のようにとります。
これからは、下の行列をRで表すことにします。
| 1 0・・・・・ |
| 0 1 ・・・・ |
| 0 0 cosθ・・・sinθ・・ | ← p行
| 0 0 0 1・・・ |
| ・・・ |
| ・・-sinθ・・cosθ・・ | ← q行
| ・・・・・・・・・・1 0 |
| ・・・・・・・・・・・1 |
↑ ↑
p列 q列
複雑に見えますがそんなことはありません。
p,q行とはこれから説明しますが、「ピボット」と言われる物です。
ピボットの部分には直交行列の基本形が入ります。
残りの部分は単位行列と同じように、
対角成分が１，残りはゼロの成分になっています。
さて、これを使って対角化するにはピボットの選び方が
問題になってくるわけですが、いくつかやり方があります。
・絶対値の大きい非対角成分を順に消す
・０でない非対角成分を適当に消していく
・絶対値がある値以上の非対角成分を適当に消し、徐々に０に近づける
ここでは、一番最初の方法を適用していきます。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★Rによる直交変換★
ここでは、「絶対値の大きい非対角成分を順に消す」という
アプローチを取ることにしました。
よって、まずはAの非対角要素の中から最大値a_pqを探し出します。
このときのa_pqがAのp行目q列目という値によって
ピボットを決定します。
こうすることによって、絶対値が大きい成分から
消していくことができます。
このようにRを作り、
A(k+1) = t[R](k)A(k)[R](k)の演算を行うことによって
固有値に収束していきます。
また、同時に単位行列UとRの積を求めることによって
固有ベクトルを求めることができます。
これは、
U(k) = U(k-1)[R](k)と書けます。
もちろんU(0)=Eですが、直交変換することによって
単位行列ではなくなります。
お気づきかと思いますが、これは掃き出し法の過程に相当するわけです。
単位行列が登場するところから想像に難くないですよね。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★実際の計算★
では、実際にt[R]A[R]の演算をしたいと思います。
ここでは便宜上、B = t[R]A[R]と書くことにします。
ここで、この演算をすると、pとqの行と列に関係する成分に変化が生じます。
実際に行列の演算を行うと、
b_pp = a_pp(cosθ)^2 + a_qq(sinθ)^2 - 2a_pq(sinθcosθ)
b_qq = a_pp(sinθ)^2 + a_qq(cosθ)^2 + 2a_pq(sinθcosθ)
b_pq = (a_pp - a_qq)(sinθcosθ) + a_pq((cosθ)^2 - (sinθ)^2)
= 1/2 (a_pp - a_qq)sin2θ + (a_pq cos2θ)
b_qp = b_pq
b_pj = a_pj cosθ - a_qj sinθ
b_qj = a_pj sinθ + a_qj cosθ
b_ip = a_ip cosθ - a_iq sinθ
b_iq = a_ip sinθ + a_iq cosθ 但しi,jはカウンタ変数
ここまでは、今まで学習してきたように
行列の積を求め、展開しただけです。
イメージの沸かない方は実際にやってみましょう。
＃メルマガでは書ききれないんで省略させてください ^^;
次に、整理することを考えます。
θは勝手に決めることができるので、b_pqをゼロにするように
θを決めることにすます。
すると、 1/2 (a_pp - a_qq)sin2θ + (a_pq cos2θ) = 0 より
∴ sinθ = √[(1-γ)/2]・(αβの符号）
cosθ = √[1-(sinθ)^2]
∵α = -a_pq
β = (a_pp - a_qq)/2
γ = |β|√[α^2 + β^2]
とまとめることができます。
また、固有ベクトルは以下のように書けます。
U_ip = U_ip cosθ - u_iq sinθ
U_iq = U_ip sinθ + u_iq cosθ
＃漸化式です。右辺に代入し、新たなUが求まります。
sinθ,cosθに関しては上で求めたものを利用できます。
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★アルゴリズムのまとめ★
[1] 対称行列Aの最大値を探し、ピボットp,qを決める
↓
[2] t[R]A[R]を求め、固有値を収束させる
U[R]を求め、固有ベクトルに収束させる
↓
[3] 最大値がゼロになる（許容相対誤差EPS以下になる）まで
手順[1][2]を繰り返す
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★サンプルプログラム★
今回学習したアルゴリズムを用い、ヤコビ法を実現します。
#include<iostream>
#include<cmath>
using namespace std;
const double EPS = 0.0001; /* 許容相対誤差 */
static int dim; /* 次数 */
double getMax(double** a,int& p,int& q){ /* 最大値求め位置p,qを返す */
double max = 0.0;
int i,j;
for(i=0; i<dim-1; i++){ /* 対角成分を検索するのでdim-1は無視 */
for(j=i+1;j<dim;j++){
if(fabs(a[i][j]) > max){
p=i;
q=j;
max=fabs(a[i][j]);
}
}
}
return max;
}
void transform(double** a,double** u,int p,int q){ /* 直交変換 */
double alpha,beta,gamma; /* α、β、γ */
double sine,cosine,w; /* sinθ,cosθ,一時値 */
double wa,wb,wc; /* 一時的にa_pp,a_pq,a_qqを格納 */
int i,j;
wa = a[p][p]; /* 前節の導出に基づき */
wb = a[p][q]; /* 必要な値を求める */
wc = a[q][q];
alpha = -wb;
beta = (wa-wc)/2;
gamma = fabs(beta)/sqrt(alpha*alpha+beta*beta);
sine = sqrt((1.0-gamma)/2);
if(alpha*beta < 0) sine = -sine;
cosine = sqrt(1.0-sine*sine);
for(j=0; j<dim; j++){
w = a[p][j]*cosine - a[q][j]*sine; /*漸化式なので分けないとNG*/
a[q][j] = a[p][j]*sine + a[q][j]*cosine;
a[p][j] = w;
}
for(j=0; j<dim; j++){
a[j][p] = a[p][j];
a[j][q] = a[q][j];
}
w = 2.0*wb*sine*cosine;
a[p][p] = wa*cosine*cosine + wc*sine*sine -w;
a[q][q] = wa*sine*sine + wc*cosine*cosine +w;
a[p][q] = a[q][p] = 0;
/* 固有ベクトルの導出 */
for(i=0; i<dim; i++){
w = u[i][p]*cosine - u[i][q]*sine; /*漸化式なので分けないとNG*/
u[i][q] = u[i][p]*sine + u[i][q]*cosine;
u[i][p] = w;
}
}
void disp(double** a,double** u){
int i,j;
for(i=0;i<dim;i++) cout << a[i][i] << ' ' << endl;
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
cout.width(10);
cout.precision(6);
cout << u[i][j] << ' ';
}
cout << endl;
}
}
int main(){
int i, j; /* カウンタ */
int p,q; /* ピボット */
double** a; /* 連立方程式の各係数 */
cout << "行列の次数は : "; /* 次元数の入力 */
cin >> dim;
a = new double*[dim]; /* メモリ確保 */
for(i=0;i<dim;i++)a[i] = new double[dim*2];
for(i=0;i<dim;i++){ /* 係数入力部 */
for(j=i;j<dim;j++){ /* 対称行列を入力 */
cout << "a[" << i << "][" << j << "] : ";
cin >> a[i][j];
a[j][i] = a[i][j];
}
}
double** u = new double*[dim]; /* Uをメモリ確保 */
for(i=0;i<dim;i++) u[i] = new double[dim]; /* 単位行列を代入 */
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
if(i==j) u[i][j] = 1;
else u[i][j] = 0;
}
}
while(1){
double max;
max = getMax(a,p,q); /* p,qにピボットが渡される */
transform(a,u,p,q); /* 直交変換 */
if(max<EPS) break; /* 収束し、終了 */
}
disp(a,u); /* 表示 */
return 0;
}
─[解説]────────────────────────────────
getMax()はAの最大値を求めます。
transform()はAとUを直交変換します。
disp()は結果の表示を行います。
また、入力は対称行列でないとまずいので、入力部に工夫があります。
内容に関しては、全て解説済みなので省略します。
コメントを良く読めば理解して頂けると思います。
────────────────────────────────────
─[実行結果]──────────────────────────────
| 3 2 -2 |
| 2 3 -2 | の固有値・固有ベクトルを求めてみましょう。
| -2 -2 8 |
行列の次数は : 3
a[0][0] : 3
a[0][1] : 2
a[0][2] : -2
a[1][1] : 3
a[1][2] : 2
a[2][2] : 8
1.05384e-16
5
9
0.666667 0.707107 -0.235702
-0.666667 0.707107 0.235702
0.333333 0 0.942809
ここから分かるのは、固有値0,5,9.
＃1e-16は0.0000000000000001のことですから、ほぼゼロですね。
固有ベクトル( 2 ) ( 1 ) ( 1 ) とわかりますね。
( -2 ) , ( 1 ) , ( 1 )
( 1 ) ( 0 ) ( 4 )
＃固有ベクトルは、その"比"に意味がありましたよね。
他にもいろいろ試してみてください。
驚くほど（？）単純に固有値・固有ベクトルを求めることができました。
＃手計算よりは遙かに楽だとおもいます。
────────────────────────────────────
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★今日のポイント★
・対称行列の固有ベクトルを縦に並べたものは、直交行列になる
・ヤコビ法は、元の行列に適当な直交行列を掛け合わせ、
直交変換することによって対角行列に近づくことを利用する
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
★予告★
複素関数について学習します
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
講読解除はこちら
http://web1.freecom.ne.jp/~mu-home/mmg/cpp.html
バックナンバーはこちら
http://web1.freecom.ne.jp/~mu-home/mmg/cpp.html
内容について質問やご意見など
smukun.com
筆者のホームページ（むーくんの理学的なんでも講座）
http://www.hello.sh/nandemo/
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊